Preliminary tests for dwarf galaxies orbits' plots with plotly
(given their ICRS coordinates)
In [1]:
from astropy import units as u
from astropy.coordinates import SkyCoord
from astropy.table import QTable
import astropy as a0
import numpy as np
import matplotlib.pyplot as plt
import math
import plotly
import plotly.graph_objs as go
import plotly.express as px
from visorbits import rotations as rot
from visorbits import ellipse as consec
plotly.offline.init_notebook_mode()
1. Initialize coordsiand transform to Galactocentric Cartesian system
In [2]:
c = {'Draco': {'ra': 15.*(17+20/60.+ 12.4/3600.)*u.deg, 'dec': (57+54/60.+55/3600.)*u.deg,\
'D': 75.8*u.kpc, 'v_r': -291.0*u.km/u.s,\
'mu_ra_cosd': 44.*10**(-3)*u.mas/u.yr, 'mu_dec': -188.*10**(-3)*u.mas/u.yr},\
'Sculptor': {'ra': 15.*(1+ 0/60.+ 9.4/3600.)*u.deg, 'dec': (-33-42/60.-33/3600.)*u.deg,\
'D': 83.9*u.kpc, 'v_r': 111.4*u.km/u.s,\
'mu_ra_cosd': 100.*10**(-3)*u.mas/u.yr, 'mu_dec': -158.*10**(-3)*u.mas/u.yr},\
'Fornax': {'ra': 15.*(2+ 39/60.+ 59.3/3600.)*u.deg, 'dec': (-34-26/60.-57/3600.)*u.deg, \
'D': 147.2*u.kpc, 'v_r': 55.3*u.km/u.s,\
'mu_ra_cosd': 381.*10**(-3)*u.mas/u.yr, 'mu_dec': -359.*10**(-3)*u.mas/u.yr},\
'LeoII': {'ra': 15.*(11+13./60.+28.8/3600.)*u.deg, 'dec': (22+ 9/60.+ 6/3600.)*u.deg, \
'D': 233.0*u.kpc, 'v_r': 78.0*u.km/u.s,\
'mu_ra_cosd': -109.*10**(-3)*u.mas/u.yr, 'mu_dec': -150.*10**(-3)*u.mas/u.yr},\
'LeoI': {'ra': 15.*(10+ 8./60.+28.1/3600.)*u.deg, 'dec': (12+18/60.+23/3600.)*u.deg,\
'D': 258.2*u.kpc, 'v_r': 282.5*u.km/u.s,\
'mu_ra_cosd': -50.*10**(-3)*u.mas/u.yr, 'mu_dec': -120.*10**(-3)*u.mas/u.yr}}
set_of_icrs = {}
for dSph in c:
set_of_icrs[dSph] = SkyCoord(ra=[c[dSph]['ra']], dec=[c[dSph]['dec']], distance=[c[dSph]['D']],\
pm_ra_cosdec=[c[dSph]['mu_ra_cosd']], pm_dec = [c[dSph]['mu_dec']],\
radial_velocity=[c[dSph]['v_r']], frame='icrs')
sGC0 = set_of_icrs[dSph].transform_to(a0.coordinates.Galactocentric)
if dSph=='Draco':
sGC = sGC0.to_table()
sGC = a0.table.hstack([sGC, QTable({'Name': [dSph]}, names=['Name'])])
else:
sGC = a0.table.vstack([sGC, a0.table.hstack([sGC0.to_table(), QTable({'Name': [dSph]}, names=['Name'])])], join_type='inner')
In [3]:
data = go.Scatter3d(x=sGC['x'], y=sGC['y'], z=sGC['z'], text = sGC['Name'], mode='markers', marker=dict(size=4), name='Dwarves')
fig = go.Figure(data)
for i, n in enumerate(sGC['Name']):
fig.add_trace(
go.Scatter3d(
x=[sGC['x'][i].value, sGC['x'][i].value+sGC['v_x'][i].value],
y=[sGC['y'][i].value, sGC['y'][i].value+sGC['v_y'][i].value],
z=[sGC['z'][i].value, sGC['z'][i].value+sGC['v_z'][i].value],
mode='lines', legendgroup=n, hovertext=n,
name="{}".format(n) ) )
fig.update_layout(title='Figure 1.',
autosize=True,
scene_camera_eye=dict(x=1.2, y=1.2, z=1.2),
scene = dict(
xaxis_title='X', yaxis_title='Y', zaxis_title='Z',
xaxis = dict(nticks=15,range=[-770,770],),
yaxis = dict(nticks=15,range=[-770,770],),
zaxis = dict(nticks=15,range=[-770,770],),),
width=550, height=550, template='none', scene_aspectmode='cube')
fig.show()
2. Prepare plotly local functions
| Function | Arguments | Output | Description |
|---|---|---|---|
| psurface |
|
|
Define our plane which contains the zero point |
| rotatedset |
|
|
Step-bystep rotation of initil dataset. |
| initfig_with_surface |
|
fig object of plotly |
Initilizes fig |
| plotvectors |
|
-- | Updates fig |
| plottraces |
fig, traces, sys_legends, colors
|
-- | Updates fig |
In [4]:
def psurface(p0, p1):
# p0 and p1 together with (0,0,0) define the plane.
# Find the normal vector:
norm = rot.vect_m(p0, p1)
norm /= np.sqrt(np.sum(norm**2))
# vtau_test is the direction of motion
# p1 marks the location of our point (which lies on the future trace)
vtau_test = rot.vect_m(norm, p1)
#tau_ort = rot.vect_m(norm_ort, r_ort) # vtau_test
tau = vtau_test/np.sqrt(np.sum(vtau_test**2))
# Definition of the final XY plane:
x = np.array([0, p0[0], p1[0]])
y = np.array([0, p0[1], p1[1]])
z0 = lambda xl, yl: (xl*norm[0] +yl*norm[1])*(-1)/norm[2]
xv, yv = np.meshgrid(x,y)
zv = z0(xv, yv)
return norm, tau, x, y, zv
def rotdataset(initial, p1, tau_ort, phi4):
# Create roational matrices
r1, r2, r3, r4 = rot.reverse_tr(p1[0], p1[1], p1[2], tau_ort, phi4=phi4)
# Application of rotational matrices, step-by-step:
ell1 = np.dot(r1, initial.T).T
ell2 = np.dot(r1, np.dot(r2, initial.T)).T
ell3 = np.dot(r1, np.dot(r2, np.dot(r3, initial.T))).T
ell4 = np.dot(np.dot(r1, np.dot(r2, np.dot(r3,r4))), initial.T).T
return initial, ell1, ell2, ell3, ell4, r1, r2, r3, r4
def initifig_with_surface(x, y, zv):
# Figure: vectors and 3-axes of coordinates.
fig = go.Figure(data=[go.Surface(x=x, y=y, z=zv,\
surfacecolor=np.full(zv.shape,'green'),opacity=0.2, showscale = False, name='Destination')])
return fig
def plotvectors(fig, vectors, vc_legends, vc_colors):
for i, v in enumerate(vectors):
fig.add_trace(go.Scatter3d(x=[v[0][0], (v[1][0]-v[0][0])*3],\
y=[v[0][1], (v[1][1]-v[0][1])*3],\
z=[v[0][2], (v[1][2]-v[0][2])*3], mode='lines',\
line = dict(color=vc_colors[i]) ,line_width=3,showlegend=True,name=vc_legends[i],))
def plottraces(fig, traces, sys_legends, colors):
for i, d in enumerate(traces):
fig.add_trace(go.Scatter3d(x=d[:,0], y=d[:,1], z=d[:,2],\
line = dict(color=colors[i]),\
marker_size=1,showlegend=True,name=sys_legends[i],))
3. Tests.
3.1 Three-axes:
In [5]:
# Initialize two points, which with (0,0,0) determine the plane of final XY
p0x = 0.201; p0y = -0.8;
p1x = -0.1; p1y = 0.85
p0 = np.array([p0x, p0y, np.sqrt(1-p0x**2-p0y**2)])
p1 = np.array([p1x, p1y, np.sqrt(1-p1x**2-p1y**2)])
initial = np.array([[0,1,0,0,0,0],\
[0,0,0,1,0,0],\
[0,0,0,0,0,1]]).T
norm, tau, x, y, zv = psurface(p0, p1)
initial, ell1, ell2, ell3, ell4, r1, r2, r3, r4 = rotdataset(initial, p1, tau, phi4=0)
fig0 = initifig_with_surface(x, y, zv)
# Prepare settings for vectors' plot:
vectors = np.array([[np.zeros(3,), norm], [np.zeros(3,), p1], [p1, p1+tau]])
vc_legends = ['Normal', 'Rad-vector of the point which belongs to the final ellipse (x3)', 'Tangential velocity direction']
vc_colors = ['pink', 'orange', 'cyan']
# Prepare settings for lines' plot:
traces = list([initial, ell1, ell2, ell3, ell4])
sys_legends = ['Initial', 'Aligned with XY-projection of point, which belongs to final',\
'Final X is he point radius-vector', 'XY is in the correct plane', 'Correct phase']
colors = ['black', 'red', 'orange', 'lime', 'cyan']
plotvectors(fig0, vectors, vc_legends, vc_colors)
plottraces(fig0, traces, sys_legends, colors)
fig0.update_layout(title='Figure 2. Test rotation of 3-axes.', autosize=True,
scene_camera_eye=dict(x=1.2, y=1.2, z=1.2),
scene = dict(xaxis_title='X', yaxis_title='Y', zaxis_title='Z',\
xaxis = dict(nticks=5), yaxis = dict(nticks=5), zaxis = dict(nticks=5),),\
width=1000, height=700, template='none', scene_aspectmode='data')
fig0.show()
3.2 Simple ellipse:
In [6]:
# Initialize two points, which with (0,0,0) determine the plane of final XY
p0x = 0.201; p0y = -0.8;
p1x = -0.1; p1y = 0.85
p0 = np.array([p0x, p0y, np.sqrt(1-p0x**2-p0y**2)])
p1 = np.array([p1x, p1y, np.sqrt(1-p1x**2-p1y**2)])
norm, tau, x, y, zv = psurface(p0, p1)
# Definition of the orbit in its own space:
# pv0 - is our curent point with phi4 phase angle from pericenter
ang = np.linspace(-np.pi, np.pi, 50);
e = 0.6
r = np.sqrt(np.sum(p1**2))/(1-e)
xcircle = r*(1-e**2)/(1+e*np.cos(ang))*np.cos(ang)
ycircle = r*(1-e**2)/(1+e*np.cos(ang))*np.sin(ang)
zcircle = np.zeros(len(ang))
points1 = np.hstack([np.reshape(xcircle, (len(xcircle), 1)), np.reshape(ycircle, (len(xcircle), 1))])
initial = np.hstack([points1, np.reshape(zcircle, (len(xcircle), 1))])
initial, ell1, ell2, ell3, ell4, r1, r2, r3, r4 = rotdataset(initial, p1, tau, phi4=0)
fig1 = initifig_with_surface(x, y, zv)
# Prepare settings for vectors' plot:
vectors = np.array([[np.zeros(3,), norm], [np.zeros(3,), p1], [p1, p1+tau]])
vc_legends = ['Normal', 'Rad-vector of the point which belongs to the final ellipse (x3)', 'Tangential velocity direction']
vc_colors = ['pink', 'orange', 'cyan']
# Prepare settings for lines' plot:
traces = list([initial, ell1, ell2, ell3, ell4])
sys_legends = ['Initial', 'Aligned with XY-projection of point, which belongs to final',\
'Final X is he point radius-vector', 'XY is in the correct plane', 'Correct phase']
colors = ['black', 'red', 'orange', 'lime', 'cyan']
plotvectors(fig1, vectors, vc_legends, vc_colors)
plottraces(fig1, traces, sys_legends, colors)
fig1.update_layout(title='Figure 3. Test with ellipses.', autosize=True,
scene_camera_eye=dict(x=1.2, y=1.2, z=1.2),
scene = dict(xaxis_title='X', yaxis_title='Y', zaxis_title='Z',\
xaxis = dict(nticks=5), yaxis = dict(nticks=5), zaxis = dict(nticks=5),),\
width=1000, height=700, template='none', scene_aspectmode='data')
fig1.show()
3.3 All rotations of our data:
In [7]:
s = sGC
s = sGC
ellipses = {}
for key in list(sGC['Name']):
ellipses[key] = {'x': [], 'y': [], 'z': [], 'set': []}
for i00, n in enumerate(s['Name']):
r_dwarf = np.array([s['x'][i00].value, s['y'][i00].value, s['z'][i00].value])
v_dwarf = np.array([s['v_x'][i00].value, s['v_y'][i00].value, s['v_z'][i00].value])
r0 = np.sum(r_dwarf**2)**0.5;
v0 = np.sum(v_dwarf**2)**0.5
r_ort = r_dwarf/r0
#norm_ort = rot.vect_m(r_ort, v_dwarf)/np.sum(rot.vect_m(r_ort, v_dwarf)**2)**0.5
norm_ort, tau_ort_wrong, x, y, zv = psurface(r_dwarf, v_dwarf)
tau_ort = rot.vect_m(norm_ort, r_ort) # vtau_test
vr = v0*np.sum(r_ort*v_dwarf)/(np.sum(r_ort**2)**0.5*np.sum(v_dwarf**2)**0.5) * r_ort
vt = tau_ort * np.sqrt(v0**2 - np.sum(vr**2))
vt0 = np.sqrt(np.sum(vt**2))
e = consec.ecc(r0*u.kpc, v0*u.km/u.s, vt0*u.km/u.s)
a = consec.a_ax(r0*u.kpc, v0*u.km/u.s)
psi_current = consec.angl_ecc(r0, e, a)
while abs(psi_current)>2*np.pi:
psi_current +=2*np.pi*(-1*psi_current)/(abs(psi_current))
phi4 = consec.angl_c(psi_current, e)
while abs(phi4)>2*np.pi:
phi4 +=2*np.pi*(-1*phi4)/(abs(phi4))
# Definition of the orbit in its own space:
# r_dwarf - is our curent point with phi4 phase angle from pericenter
ang = np.linspace(0, np.pi*2, 50);
xcircle = a*(1-e**2)/(1+e*np.cos(ang))*np.cos(ang)
ycircle = a*(1-e**2)/(1+e*np.cos(ang))*np.sin(ang)
zcircle = np.zeros(len(ang))
points1 = np.hstack([np.reshape(xcircle, (len(xcircle), 1)), np.reshape(ycircle, (len(xcircle), 1))])
initial = np.hstack([points1, np.reshape(zcircle, (len(xcircle), 1))])
p1 = np.array([xcircle[0], ycircle[0], zcircle[0]])
p2 = np.array([a*(1-e**2)/(1+e*np.cos(phi4))*np.cos(phi4),\
a*(1-e**2)/(1+e*np.cos(phi4))*np.sin(phi4), 0.])
# Create roational matrices
initial, ell1, ell2, ell3, ell4, r1, r2, r3, r4 = rotdataset(initial, r_dwarf, vt, phi4=phi4)
ellipses[n]['x'] = ell4[:,0]; ellipses[n]['y'] = ell4[:,1]; ellipses[n]['z'] = ell4[:,2]; ellipses[n]['set'] = ell4
p1 = np.dot(np.dot(r1, np.dot(r2, np.dot(r3,r4))), p1)
p2 = np.dot(np.dot(r1, np.dot(r2, np.dot(r3,r4))), p2)
fig2 = go.Figure()
# Prepare settings for vectors' plot:
vectors = np.array([[np.zeros(3,), norm_ort*a], [np.zeros(3,), r_dwarf], [r_dwarf, r_dwarf+v_dwarf],\
[np.zeros(3,), p1], [np.zeros(3,), p2]])
vc_legends = ['Normal', 'Rad-vector of the point which belongs to the final ellipse (x3)',\
f'Velocity direction, phase angle is {phi4*180./np.pi}', 'Pericenter', 'The Dwarf']
vc_colors = ['pink', 'orange', 'cyan', 'black', 'red']
# Prepare settings for lines' plot:
traces = list([initial, ell1, ell2, ell3, ell4])
ell_legends = ['Initial', 'Aligned with XY-projection of point, which belongs to final',\
'Final X is he point radius-vector', 'XY is in the correct plane', 'Correct phase: X is in the pericenter']
colors = ['black', 'red', 'orange', 'lime', 'cyan']
plotvectors(fig2, vectors, vc_legends, vc_colors)
plottraces(fig2, traces, ell_legends, colors)
fig2.update_layout(title=f'Figure {4+i00}. {n}', autosize=True,
scene_camera_eye=dict(x=1.2, y=1.2, z=1.2),
scene = dict(xaxis_title='X', yaxis_title='Y', zaxis_title='Z',\
xaxis = dict(nticks=5), yaxis = dict(nticks=5), zaxis = dict(nticks=5),),\
width=1000, height=700, template='none', scene_aspectmode='data')
fig2.show()
3.4 Plot everything:
In [8]:
fig3 = go.Figure()
colors_dwarf = ['red', 'orange', 'green', 'blue', 'purple']
for i, n in enumerate(list(sGC['Name'])):
d = ellipses[n]['set']
plottraces(fig3, list([d]), [n], [colors_dwarf[i]])
fig3.add_trace(
go.Scatter3d(
x=[sGC['x'][i].value, sGC['x'][i].value+sGC['v_x'][i].value],
y=[sGC['y'][i].value, sGC['y'][i].value+sGC['v_y'][i].value],
z=[sGC['z'][i].value, sGC['z'][i].value+sGC['v_z'][i].value],
mode='lines', legendgroup=n, hovertext=n, line=dict(color=colors_dwarf[i]),
name="{} position".format(n) ) )
fig3.update_layout(title='Figure..', autosize=True,
scene_camera_eye=dict(x=1.2, y=1.2, z=1.2),
scene = dict(xaxis_title='X', yaxis_title='Y', zaxis_title='Z',\
xaxis = dict(nticks=5), yaxis = dict(nticks=5), zaxis = dict(nticks=5),),\
width=1000, height=700, template='none', scene_aspectmode='data')
fig3.show()
In [ ]: